Ce projet s’inscrit dans le cadre de l’étude de la planète Mercure, en particulier de sa géochimie de surface.
Grâce aux données de la sonde MESSENGER et à des fichiers satellites (BMP, DAT, TIF, CSV), l’objectif est de construire un tableau 3D qui rassemble différentes couches géochimiques (Mg/Si, Al/Si, Fe/Si, etc.).

Le but final est de produire des cartes géochimiques de la surface de Mercure, afin de comparer la distribution des éléments et mieux comprendre l’évolution magmatique de la planète.



1 Préparation des données

Avant de pouvoir analyser et cartographier les données, il faut :
1. Nettoyer l’environnement R pour éviter les conflits,
2. Installer et charger les packages nécessaires (visualisation, traitement raster, tables interactives…),
3. Définir les chemins vers les dossiers où se trouvent les données,
4. Vérifier l’existence du dossier de données pour garantir la reproductibilité.

Cette étape assure une base solide pour la suite du traitement.

1.1 packages & chemins

Ici, on installe automatiquement les packages manquants, puis on les charge silencieusement.
On définit ensuite le dossier de données (mercury1440x720) et le fichier de sortie (result_array.rds).
Enfin, on installe et charge le package mercure.ulg, développé pour traiter ce type de données.

# Nettoyage minimal (optionnel)
# rm(list = ls())
options(repos = c(CRAN = "https://cloud.r-project.org"))
# Packages requis (on installe si manquants)
need <- c("plotly", "raster", "bmp", "devtools", "rstudioapi", "DT", "viridisLite")
to_install <- setdiff(need, rownames(installed.packages()))
if (length(to_install)) install.packages(to_install, quiet = TRUE)

# Chargement
invisible(lapply(need, require, character.only = TRUE))

# Base path robuste (RStudio ou knitr)
base_path <- tryCatch(
  dirname(rstudioapi::getSourceEditorContext()$path),
  error = function(e) getwd()
)

# Dossier de données (à adapter si besoin)
directory_path <- file.path(base_path, "mercury1440x720")
result_file    <- file.path(base_path, "result_array.rds")

# Installation/chargement du package d'analyse (GitHub)
if (!requireNamespace("mercure.ulg", quietly = TRUE)) {
  devtools::install_github("Zaeron1/mercure.ulg", upgrade = "never", quiet = TRUE)
}
library(mercure.ulg)

# Sécurité : le dossier de données doit exister
if (!dir.exists(directory_path)) {
  stop(paste("Le dossier de données est introuvable :", directory_path))
}

1.2 Construction du tableau 3D

Pour analyser conjointement toutes les données satellites, il est nécessaire de les regrouper dans une seule structure : un tableau 3D.
Chaque couche du tableau correspond à un fichier de données (par exemple une carte Mg/Si ou un modèle d’élévation).

La fonction process_files_to_3d_array() :
- parcourt tous les fichiers du dossier,
- applique les facteurs de correction nécessaires pour normaliser les valeurs,
- redimensionne les images si besoin,
- renvoie un tableau 3D où chaque plan correspond à une carte.

Cette fonction est donc au cœur du projet, car elle permet d’unifier tous les formats (BMP, DAT, TIF, CSV) dans une seule base exploitable.

process_files_to_3d_array <- function(directory_path, correction_factors, target_dim = c(720, 1440)) {
  file_list   <- list.files(directory_path, full.names = TRUE)
  n_files     <- length(file_list)
  if (n_files == 0) stop("Aucun fichier détecté dans: ", directory_path)

  array_3d    <- array(0, dim = c(target_dim[1], target_dim[2], n_files))
  layer_names <- basename(file_list)

  apply_correction <- function(mat, factor) {
    if (is.null(factor)) factor <- 1
    mat * factor
  }

  for (i in seq_along(file_list)) {
    file_path <- file_list[i]
    file_name <- layer_names[i]
    ext <- tolower(tools::file_ext(file_name))

    if (ext == "bmp") {
      bmp_matrix <- bmp::read.bmp(file_path)
      array_3d[,,i] <- apply_correction(bmp_matrix, correction_factors[[file_name]])

    } else if (ext == "dat") {
      dat <- as.matrix(read.table(file_path, header = FALSE))
      array_3d[,,i] <- mercure.ulg::from03602180(dat)

    } else if (ext == "tif") {
      tif <- raster::raster(file_path)
      fact <- c(max(1, ncol(tif) / target_dim[2]),
                max(1, nrow(tif) / target_dim[1]))
      tif_small <- raster::aggregate(tif, fact = fact)
      mat <- as.matrix(tif_small, ncol = target_dim[2], nrow = target_dim[1])
      array_3d[,,i] <- mat * 0.5  # correction simple d'exemple

    } else if (ext == "csv") {
      lines <- gsub(",\\s+", ",", gsub("(\\d)\\s+(\\d)", "\\1,\\2", readLines(file_path)))
      tmp <- tempfile(fileext = ".csv")
      writeLines(lines, tmp)
      donnees <- read.csv(tmp, header = FALSE)
      unlink(tmp)
      array_3d[,,i] <- mercure.ulg::from03602180(as.matrix(donnees[,-1]))

    } else {
      warning("Type non supporté: ", ext, " pour le fichier ", file_name)
    }
  }

  # Info couches (utile pour navigation)
  layer_info <- data.frame(Index = seq_along(layer_names), Layer = layer_names)
  attr(array_3d, "layer_info") <- layer_info
  array_3d
}

1.3 Construction/exécution + résumé visuel

Une fois la fonction définie, on peut l’exécuter avec les fichiers satellites de Mercure.
Cette étape se décompose en plusieurs sous-parties :
1. Définition et affichage des facteurs de correction,
2. Construction du tableau 3D,
3. Résumés rapides (dimensions et taille mémoire),
4. Génération de la liste des couches disponibles.


1.3.1 🔹 Facteurs de correction

Les fichiers BMP bruts sont normalisés entre 0 et 255.
Afin de convertir ces valeurs en rapports géochimiques (ex. Mg/Si, Al/Si), on applique un facteur de correction spécifique à chaque fichier.
Ces facteurs sont listés dans le tableau ci-dessous :

correction_factors <- list(
  "mgsi.bmp" = 0.860023 / 255.0,
  "alsi.bmp" = 0.402477 / 255.0,
  "ssi.bmp"  = 0.161680 / 255.0,
  "fesi.bmp" = 0.117737 / 255.0,
  "casi.bmp" = 0.318000 / 255.0
)

# Transformation en tableau
correction_table <- data.frame(
  Fichier = names(correction_factors),
  Facteur = unlist(correction_factors)
)

# Affichage lisible (interactif si DT dispo)
if (requireNamespace("DT", quietly = TRUE)) {
  DT::datatable(correction_table,
                rownames = FALSE,
                options = list(dom = 't'),
                caption = "Facteurs de correction appliqués aux fichiers BMP")
} else {
  knitr::kable(correction_table, digits = 6,
               caption = "Facteurs de correction appliqués aux fichiers BMP")
}

1.3.2 🔹 Construction du tableau 3D

La fonction process_files_to_3d_array() est appliquée au dossier mercury1440x720. Chaque fichier est lu, corrigé et stocké dans un tableau 3D, où : • la dimension 1 correspond aux lignes (latitude), • la dimension 2 correspond aux colonnes (longitude), • la dimension 3 correspond aux différentes couches géochimiques ou géophysiques.

result_array <- process_files_to_3d_array(directory_path, correction_factors)
saveRDS(result_array, file = result_file) # Sauvegarder le tableau 3D dans un fichier RDS
#result_array <- readRDS(result_file) # Lire le tableau 3D depuis le fichier RDS
rm(list=ls()[!ls() %in% c("result_array")]) # Supprimer toutes les variables sauf le tableau 3D

1.3.3 Liste de couches disponibles

Chaque couche correspond à un fichier source (BMP, DAT, TIF, CSV). Le tableau ci-dessous récapitule les couches disponibles dans result_array.

layer_info <- attr(result_array, "layer_info")

if (requireNamespace("DT", quietly = TRUE)) {
  DT::datatable(layer_info,
                rownames = FALSE,
                options = list(pageLength = 10),
                caption = "Liste des couches présentes dans le tableau 3D")
} else {
  knitr::kable(layer_info,
               caption = "Liste des couches présentes dans le tableau 3D")
}

1.3.3.1 Obtention d’une couche spécifique

Pour accéder à une couche particulière, on peut utiliser son index ou son nom de fichier.

layer_matrix <- get_layer_as_matrix(result_array, 2) # EXEMPLE: Extraire une couche spécifique (index 2 ici) du tableau 3D sous forme de matrice

IMPORTANT: le tableau généré dans la console renseigne les index qu’il faut utiliser dans les prochaines fonctions poru extraire les bonnes matrices du 3D array

2 Résultats des analyses

Avec le tableau 3D construit, on peut désormais extraire des couches spécifiques pour analyse et visualisation.
Par exemple, on peut extraire la couche DEM (index 4) et la visualiser avec une palette de couleurs adaptée

2.1 Digital Elevation Model (DEM) de Mercure

plot_matrix <- function(layer_index_or_name, main_title, legend_name, color_palette) {
  # check si l'argument est une matrice dans l'environment ou un index à extraire de result_array
  if (is.matrix(layer_index_or_name)) {
    layer_matrix <- layer_index_or_name
  } else if (is.numeric(layer_index_or_name)) {
    if (layer_index_or_name < 1 || layer_index_or_name > dim(result_array)[3]) {
      stop("Layer index out of bounds.")
    }
    layer_matrix <- get_layer_as_matrix(result_array, layer_index_or_name)
  } else {
    stop("The input must be either a matrix or a numeric index of a layer.")
  }

  # Create a plotly heatmap
  plot_ly(
    z = layer_matrix,
    type = "heatmap",
    colors = colorRamp(c("blue", "cyan", "yellowgreen", "yellow", "orange", "orangered", "darkred")),
    showscale = TRUE,
    colorbar = list(title = legend_name)
  ) %>%
    layout(
      title = main_title,
      xaxis = list(
        title = "Longitude",
        tickvals = seq(0, ncol(layer_matrix), by = ncol(layer_matrix) / 12),
        ticktext = paste0(seq(-180, 180, length.out = 13), "°")
      ),
      yaxis = list(
        title = "Latitude",
        tickvals = seq(0, nrow(layer_matrix), by = nrow(layer_matrix) / 12),
        ticktext = paste0(seq(90, -90, length.out = 13), "°"), autorange = "reversed"
      )
    )
}
# Affichage de la couche 4 (DEM) en heatmap
plot_matrix(
  layer_index_or_name = 4, 
  main_title = "Carte du DEM de Mercure 🌍",
  legend_name = "Élévation (m)",
  color_palette = NULL
)

2.2 Classifcation en déciles géochimiques

La classification en déciles permet de segmenter les données géochimiques en 10 catégories égales.
Cela facilite l’interprétation des variations spatiales des rapports élémentaires (ex. Mg/Si).
La fonction classify_matrix() prend en entrée une couche spécifique et renvoie une matrice classifiée.

classify_matrix <- function(layer_index){
  if (layer_index < 1 || layer_index > dim(result_array)[3]) {
    stop("Layer index out of bounds.")
  }
  layer_matrix <- get_layer_as_matrix(result_array, layer_index)
  quantile <- quantile(layer_matrix, probs = seq(0, 1, by = 0.1))  # Calcul des quantile
  labels <- c("Lowest","Very Very Low", "Very Low", "Low", "Medium Low", "Medium", "Medium High", "High", "Very High", "Very Very High", "Highest")  # Étiquettes de classification
  classes <- cut(as.vector(layer_matrix), breaks = c(quantile, Inf), labels = labels, include.lowest = TRUE)  # Classification des valeurs
  class_matrix <- matrix(as.numeric(classes), nrow = nrow(layer_matrix), byrow = FALSE)  # Conversion en matrice numérique
  return(class_matrix)
}